﻿//TODO: Canny: test

using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Drawing;
using System.Drawing.Imaging;
using System.Linq;
using System.Reflection;
using System.Windows.Forms;
using LowLevelGraphics.Extensions;

//TODO: CannyEdgeDetector, TEST! DOESN'T WORK RIGHT NOW 

namespace LowLevelGraphics.Filter
{
    /// <summary>
    /// Canny edge detector
    /// </summary>
    public class CannyEdgeDetector : BaseImageFilter
    {
        // statics
        private static float GAUSSIAN_CUT_OFF = 0.005f;
        private static float MAGNITUDE_SCALE = 100F;
        private static float MAGNITUDE_LIMIT = 1000F;
        private static int MAGNITUDE_MAX = (int)(MAGNITUDE_SCALE * MAGNITUDE_LIMIT);

        // fields
        protected int m_nHeight = -1;
        protected int m_nWidth = -1;
        protected int m_nPicsize = -1;
        protected int[] m_aData = null;
        protected int[] m_aMagnitude = null;
        protected UnsafeBitmap m_SourceImage = null;
        protected UnsafeBitmap m_EdgesImage = null;

        protected float m_fGaussianKernelRadius = -1.0f;
        protected float m_fLowThreshold = -1.0f;
        protected float m_fHighThreshold = -1.0f;
        protected int m_nGaussianKernelWidth = -1;
        protected bool m_bContrastNormalized = false;

        protected float[] m_aXConv = null;
        protected float[] m_aYConv = null;
        protected float[] m_xGradient = null;
        protected float[] m_yGradient = null;

        // constructors

        /// <summary>
        /// Constructs a new detector with default parameters.        
        /// </summary>
        public CannyEdgeDetector()
        {
            m_fLowThreshold = 2.5f;
            m_fHighThreshold = 7.5f;
            m_fGaussianKernelRadius = 2f;
            m_nGaussianKernelWidth = 16;
            m_bContrastNormalized = false;
        }

        /// <summary>
        /// used by clone to create a deep copy
        /// </summary>
        /// <param name="_cannyEdgeDetector"></param>
        internal CannyEdgeDetector(CannyEdgeDetector _cannyEdgeDetector)
        {
            this.m_fLowThreshold = _cannyEdgeDetector.m_fLowThreshold;
            this.m_fHighThreshold = _cannyEdgeDetector.m_fHighThreshold;
            this.m_fGaussianKernelRadius = _cannyEdgeDetector.m_fGaussianKernelRadius;
            this.m_bContrastNormalized = _cannyEdgeDetector.m_bContrastNormalized;
        }

        // accessors

        /// <summary>
        /// The image that provides the luminance data used by this detector to
        /// generate edges.
        /// Specifies the image that will provide the luminance data in which edges
        /// will be detected. A source image must be set before the process method
        /// is called.
        /// </summary>
        /// <returns>return the source image, or null</returns>
        public UnsafeBitmap SourceImage
        {
            get { return m_SourceImage; }
            set { m_SourceImage = value; }
        }

        /// <summary>        
        ///  Obtains an image containing the edges detected during the last call to
        ///  the process method. The buffered image is an opaque image of type
        ///  BufferedImage.TYPE_INT_ARGB in which edge pixels are white and all other
        ///  pixels are black.          
        ///  @return an image containing the detected edges, or null if the process
        ///  method has not yet been called.
        /// Sets the edges image. Calling this method will not change the operation
        /// of the edge detector in any way. It is intended to provide a means by
        /// which the memory referenced by the detector object may be reduced.
        /// </summary>
        /// <returns></returns>
        public UnsafeBitmap EdgesImage
        {
            get { return m_EdgesImage; }
            set { m_EdgesImage = value; }
        }

        /// <summary>
        /// The low threshold for hysteresis. The default value is 2.5	/// </summary>
        /// Sets the low threshold for hysteresis. Suitable values for this parameter
        /// must be determined experimentally for each application. It is nonsensical
        /// (though not prohibited) for this value to exceed the high threshold value.
        /// <returns>low hysteresis threshold</returns>
        public float LowThreshold
        {
            get { return m_fLowThreshold; }
            set
            {
                if (m_fLowThreshold < 0) throw new ArgumentException();
                m_fLowThreshold = value;
            }
        }

        /// <summary>
        /// The high threshold for hysteresis. The default value is 7.5.
        /// Sets the high threshold for hysteresis. Suitable values for this
        /// parameter must be determined experimentally for each application. It is
        /// nonsensical (though not prohibited) for this value to be less than the
        /// low threshold value.
        /// </summary>
        /// <returns>return the high hysteresis threshold</returns>
        public float HighThreshold
        {
            get { return m_fHighThreshold; }
            set
            {
                if (m_fHighThreshold < 0) throw new ArgumentException();
                m_fHighThreshold = value;
            }
        }

        /// <summary>
        /// The number of pixels across which the Gaussian kernel is applied.
        /// The default value is 16.
        /// The number of pixels across which the Gaussian kernel is applied.
        /// This implementation will reduce the radius if the contribution of pixel
        /// values is deemed negligable, so this is actually a maximum radius.
        /// </summary>
        /// <returns>return the radius of the convolution operation in pixels</returns>
        public int GaussianKernelWidth
        {
            get
            {
                return m_nGaussianKernelWidth; 
            }
            set
            {
                if (m_nGaussianKernelWidth < 2) throw new ArgumentException();
                this.m_nGaussianKernelWidth = value;
            }
        }

        /// <summary>
        /// The radius of the Gaussian convolution kernel used to smooth the source
        /// image prior to gradient calculation. The default value is 16.        
        /// Sets the radius of the Gaussian convolution kernel used to smooth the
        /// source image prior to gradient calculation.                  
        /// </summary>
        /// <returns>return the Gaussian kernel radius in pixels</returns>
        public float GaussianKernelRadius
        {
            get { return m_fGaussianKernelRadius; }
            set
            {
                if (m_fGaussianKernelRadius < 0.1f) throw new ArgumentException();
                this.m_fGaussianKernelRadius = value;
            }
        }

        /// <summary>
        /// Whether the luminance data extracted from the source image is normalized
        /// by linearizing its histogram prior to edge extraction. The default value
        /// is false.
        /// Sets whether the contrast is normalized
        /// </summary>
        /// <returns>return whether the contrast is normalized</returns>
        public bool IsContrastNormalized
        {
            get { return m_bContrastNormalized; }
            set { this.m_bContrastNormalized = value; }
        }

        public override UnsafeBitmap Execute(UnsafeBitmap _bitmap)
        {
            m_SourceImage = _bitmap;
            Process();
            m_SourceImage.Dispose();
            return _bitmap;
        }

        // methods

        /// <summary>
        /// Process will execute the real canny edge detector
        /// </summary>
        public void Process()
        {
            m_nWidth = m_SourceImage.Width;
            m_nHeight = m_SourceImage.Height;
            m_nPicsize = m_nWidth * m_nHeight;
            InitArrays();
            ReadLuminance();
            if (m_bContrastNormalized) NormalizeContrast();
            ComputeGradients(m_fGaussianKernelRadius, m_nGaussianKernelWidth);
            int nLow = (int)Math.Round(m_fLowThreshold * MAGNITUDE_SCALE);
            int nHigh = (int)Math.Round(m_fHighThreshold * MAGNITUDE_SCALE);
            PerformHysteresis(nLow, nHigh);
            ThresholdEdges();
            WriteEdges(m_aData);
        }

        // private utility methods

        /// <summary>
        /// Initialize the arrays
        /// </summary>
        private void InitArrays()
        {
            if (m_aData == null || m_nPicsize != m_aData.Length)
            {
                m_aData = new int[m_nPicsize];
                m_aMagnitude = new int[m_nPicsize];

                m_aXConv = new float[m_nPicsize];
                m_aYConv = new float[m_nPicsize];
                m_xGradient = new float[m_nPicsize];
                m_yGradient = new float[m_nPicsize];
            }
        }

        //NOTE: The elements of the method below (specifically the technique for
        //non-maximal suppression and the technique for gradient computation)
        //are derived from an implementation posted in the following forum (with the
        //clear intent of others using the code):
        //  http://forum.java.sun.com/thread.jspa?threadID=546211&start=45&tstart=0
        //My code effectively mimics the algorithm exhibited above.
        //Since I don't know the providence of the code that was posted it is a
        //possibility (though I think a very remote one) that this code violates
        //someone's intellectual property rights. If this concerns you feel free to
        //contact me for an alternative, though less efficient, implementation.

        /// <summary>
        /// Computes the gradients.
        /// </summary>
        /// <param name="kernelRadius">The kernel radius.</param>
        /// <param name="kernelWidth">Width of the kernel.</param>
        private void ComputeGradients(float kernelRadius, int kernelWidth)
        {
            //generate the gaussian convolution masks
            float[] kernel = new float[kernelWidth];
            float[] diffKernel = new float[kernelWidth];
            int kwidth = 0;
            for (kwidth = 0; kwidth < kernelWidth; kwidth++)
            {
                float g1 = Gaussian(kwidth, kernelRadius);
                if (g1 <= GAUSSIAN_CUT_OFF && kwidth >= 2) break;
                float g2 = Gaussian(kwidth - 0.5f, kernelRadius);
                float g3 = Gaussian(kwidth + 0.5f, kernelRadius);
                kernel[kwidth] = (g1 + g2 + g3) / 3f / (2f * (float)Math.PI * kernelRadius * kernelRadius);
                diffKernel[kwidth] = g3 - g2;
            }

            int initX = kwidth - 1;
            int maxX = m_nWidth - (kwidth - 1);
            int initY = m_nWidth * (kwidth - 1);
            int maxY = m_nWidth * (m_nHeight - (kwidth - 1));

            //perform convolution in x and y directions
            int x = 0;
            int y = 0;

            for (x = initX; x < maxX; x++)
            {
                for (y = initY; y < maxY; y += m_nWidth)
                {
                    int index = x + y;
                    float sumX = m_aData[index] * kernel[0];
                    float sumY = sumX;
                    int xOffset = 1;
                    int yOffset = m_nWidth;
                    for (; xOffset < kwidth ;)
                    {
                        sumY += kernel[xOffset] * (m_aData[index - yOffset] + m_aData[index + yOffset]);
                        sumX += kernel[xOffset] * (m_aData[index - xOffset] + m_aData[index + xOffset]);
                        yOffset += m_nWidth;
                        xOffset++;
                    }

                    m_aYConv[index] = sumY;
                    m_aXConv[index] = sumX;
                }
            }

            for (x = initX; x < maxX; x++)
            {
                for (y = initY; y < maxY; y += m_nWidth)
                {
                    float sum = 0f;
                    int index = x + y;
                    for (int i = 1; i < kwidth; i++)
                        sum += diffKernel[i] * (m_aYConv[index - i] - m_aYConv[index + i]);

                    m_xGradient[index] = sum;
                }
            }

            for (x = kwidth; x < m_nWidth - kwidth; x++)
            {
                for (y = initY; y < maxY; y += m_nWidth)
                {
                    float sum = 0.0f;
                    int index = x + y;
                    int yOffset = m_nWidth;
                    for (int i = 1; i < kwidth; i++)
                    {
                        sum += diffKernel[i] * (m_aXConv[index - yOffset] - m_aXConv[index + yOffset]);
                        yOffset += m_nWidth;
                    }

                    m_yGradient[index] = sum;
                }

            }

            initX = kwidth;
            maxX = m_nWidth - kwidth;
            initY = m_nWidth * kwidth;
            maxY = m_nWidth * (m_nHeight - kwidth);
            for (x = initX; x < maxX; x++)
            {
                for (y = initY; y < maxY; y += m_nWidth)
                {
                    int index = x + y;
                    int indexN = index - m_nWidth;
                    int indexS = index + m_nWidth;
                    int indexW = index - 1;
                    int indexE = index + 1;
                    int indexNW = indexN - 1;
                    int indexNE = indexN + 1;
                    int indexSW = indexS - 1;
                    int indexSE = indexS + 1;

                    float xGrad = m_xGradient[index];
                    float yGrad = m_yGradient[index];
                    float gradMag = hypot(xGrad, yGrad);

                    //perform non-maximal supression
                    float nMag = hypot(m_xGradient[indexN], m_yGradient[indexN]);
                    float sMag = hypot(m_xGradient[indexS], m_yGradient[indexS]);
                    float wMag = hypot(m_xGradient[indexW], m_yGradient[indexW]);
                    float eMag = hypot(m_xGradient[indexE], m_yGradient[indexE]);
                    float neMag = hypot(m_xGradient[indexNE], m_yGradient[indexNE]);
                    float seMag = hypot(m_xGradient[indexSE], m_yGradient[indexSE]);
                    float swMag = hypot(m_xGradient[indexSW], m_yGradient[indexSW]);
                    float nwMag = hypot(m_xGradient[indexNW], m_yGradient[indexNW]);
                    float tmp;
                    /*
                     * An explanation of what's happening here, for those who want
                     * to understand the source: This performs the "non-maximal
                     * supression" phase of the Canny edge detection in which we
                     * need to compare the gradient magnitude to that in the
                     * direction of the gradient; only if the value is a local
                     * maximum do we consider the point as an edge candidate.
                     * 
                     * We need to break the comparison into a number of different
                     * cases depending on the gradient direction so that the
                     * appropriate values can be used. To avoid computing the
                     * gradient direction, we use two simple comparisons: first we
                     * check that the partial derivatives have the same sign (1)
                     * and then we check which is larger (2). As a consequence, we
                     * have reduced the problem to one of four identical cases that
                     * each test the central gradient magnitude against the values at
                     * two points with 'identical support'; what this means is that
                     * the geometry required to accurately interpolate the magnitude
                     * of gradient function at those points has an identical
                     * geometry (upto right-angled-rotation/reflection).
                     * 
                     * When comparing the central gradient to the two interpolated
                     * values, we avoid performing any divisions by multiplying both
                     * sides of each inequality by the greater of the two partial
                     * derivatives. The common comparand is stored in a temporary
                     * variable (3) and reused in the mirror case (4).
                     * 
                     */
                    if (xGrad * yGrad <= (float)0 /*(1)*/
                        ? Math.Abs(xGrad) >= Math.Abs(yGrad) /*(2)*/
                            ? (tmp = Math.Abs(xGrad * gradMag)) >= Math.Abs(yGrad * neMag - (xGrad + yGrad) * eMag) /*(3)*/
                                && tmp > Math.Abs(yGrad * swMag - (xGrad + yGrad) * wMag) /*(4)*/
                            : (tmp = Math.Abs(yGrad * gradMag)) >= Math.Abs(xGrad * neMag - (yGrad + xGrad) * nMag) /*(3)*/
                                && tmp > Math.Abs(xGrad * swMag - (yGrad + xGrad) * sMag) /*(4)*/
                        : Math.Abs(xGrad) >= Math.Abs(yGrad) /*(2)*/
                            ? (tmp = Math.Abs(xGrad * gradMag)) >= Math.Abs(yGrad * seMag + (xGrad - yGrad) * eMag) /*(3)*/
                                && tmp > Math.Abs(yGrad * nwMag + (xGrad - yGrad) * wMag) /*(4)*/
                            : (tmp = Math.Abs(yGrad * gradMag)) >= Math.Abs(xGrad * seMag + (yGrad - xGrad) * sMag) /*(3)*/
                                && tmp > Math.Abs(xGrad * nwMag + (yGrad - xGrad) * nMag) /*(4)*/
                        )
                    {
                        m_aMagnitude[index] = gradMag >= MAGNITUDE_LIMIT ? MAGNITUDE_MAX : (int)(MAGNITUDE_SCALE * gradMag);
                        //NOTE: The orientation of the edge is not employed by this
                        //implementation. It is a simple matter to compute it at
                        //this point as: Math.atan2(yGrad, xGrad);
                    }
                    else
                    {
                        m_aMagnitude[index] = 0;
                    }
                }
            }
        }

        /// <summary>
        ///NOTE: It is quite feasible to replace the implementation of this method
        ///with one which only loosely approximates the hypot function. I've tested
        ///simple approximations such as Math.abs(x) + Math.abs(y) and they work fine.
        /// </summary>
        /// <param name="x"></param>
        /// <param name="y"></param>
        /// <returns></returns>
        private float hypot(float x, float y)
        {
            return (float)Hypot(x, y);
        }

        /// <summary>
        /// This returns the hypothenuse of a triangle
        /// </summary>
        /// <param name="_x"></param>
        /// <param name="_y"></param>
        /// <returns></returns>
        private double Hypot(double _x, double _y)
        {
            return Math.Sqrt(_x * _x + _y * _y);
        }

        /// <summary>
        /// Gaussians the specified x.
        /// </summary>
        /// <param name="x">The x.</param>
        /// <param name="sigma">The sigma.</param>
        /// <returns></returns>
        private float Gaussian(float x, float sigma)
        {
            return (float)Math.Exp(-(x * x) / (2f * sigma * sigma));
        }

        /// <summary>
        /// Performs the hysteresis.
        /// </summary>
        /// <param name="low">The low.</param>
        /// <param name="high">The high.</param>
        private void PerformHysteresis(int low, int high)
        {
            //NOTE: this implementation reuses the data array to store both
            //luminance data from the image, and edge intensity from the processing.
            //This is done for memory efficiency, other implementations may wish
            //to separate these functions.
            //Arrays.fill(data, 0);

            int nOffset = 0;
            
            int x = 0;
            int y = 0;

            for (x = 0; x < m_nWidth; x++)
            {
                for (y = 0; y < m_nHeight; y++)
                {
                    if (m_aData[nOffset] == 0 && m_aMagnitude[nOffset] >= high)
                    {
                        Follow(x, y, nOffset, low);
                    }
                    nOffset++;
                }
            }
        }

        /// <summary>
        /// Follows the specified x1.
        /// </summary>
        /// <param name="x1">The x1.</param>
        /// <param name="y1">The y1.</param>
        /// <param name="i1">The i1.</param>
        /// <param name="threshold">The threshold.</param>
        private void Follow(int x1, int y1, int i1, int threshold)
        {
            int x0 = x1 == 0 ? x1 : x1 - 1;
            int x2 = x1 == m_nWidth - 1 ? x1 : x1 + 1;
            int y0 = y1 == 0 ? y1 : y1 - 1;
            int y2 = y1 == m_nHeight - 1 ? y1 : y1 + 1;

            m_aData[i1] = m_aMagnitude[i1];
            for (int x = x0; x <= x2; x++)
            {
                for (int y = y0; y <= y2; y++)
                {
                    int i2 = x + y * m_nWidth;
                    if ((y != y1 || x != x1)
                        && m_aData[i2] == 0
                        && m_aMagnitude[i2] >= threshold)
                    {
                        Follow(x, y, i2, threshold);
                        return;
                    }
                }
            }
        }

        /// <summary>
        /// Thresholds the edges.
        /// </summary>
        private void ThresholdEdges()
        {
            for (int i = 0; i < m_nPicsize; i++)
            {
                m_aData[i] = m_aData[i] > 0 ? 0xFFFFFF : 0; //oldcode: -1 : 0xFF000000
            }
        }

        /// <summary>
        /// Luminances the specified r.
        /// </summary>
        /// <param name="r">The r.</param>
        /// <param name="g">The g.</param>
        /// <param name="b">The b.</param>
        /// <returns></returns>
        private int Luminance(float r, float g, float b)
        {
            return (int)Math.Round(0.299f * r + 0.587f * g + 0.114f * b);
        }

        /// <summary>
        /// Reads the luminance.
        /// </summary>
        private void ReadLuminance()
        {
            //int type = sourceImage.getType();
            //if (type == BufferedImage.TYPE_INT_RGB || type == BufferedImage.TYPE_INT_ARGB)
            {
                int i = 0;
                Color color = Color.Transparent;
                for (int y = 0; y < m_nHeight; y++)
                {
                    for (int x = 0; x < m_nWidth; x++)
                    {
                        int p = m_SourceImage.GetPixelDirectInt(x, y);
                        int r = (p & 0xff0000) >> 16;
                        int g = (p & 0xff00) >> 8;
                        int b = p & 0xff;
                        m_aData[i] = (int)Luminance(r, g, b); //luminance(r, g, b);
                        i++;
                    }
                }
            }
            /*else if (type == BufferedImage.TYPE_BYTE_GRAY)
            {
                byte[] pixels = (byte[])sourceImage.getData().getDataElements(0, 0, width, height, null);
                for (int i = 0; i < picsize; i++)
                {
                    data[i] = (pixels[i] & 0xff);
                }
            }
            else if (type == BufferedImage.TYPE_USHORT_GRAY)
            {
                short[] pixels = (short[])sourceImage.getData().getDataElements(0, 0, width, height, null);
                for (int i = 0; i < picsize; i++)
                {
                    data[i] = (pixels[i] & 0xffff) / 256;
                }
            }
            else if (type == BufferedImage.TYPE_3BYTE_BGR)
            {
                byte[] pixels = (byte[])sourceImage.getData().getDataElements(0, 0, width, height, null);
                int offset = 0;
                for (int i = 0; i < picsize; i++)
                {
                    int b = pixels[offset++] & 0xff;
                    int g = pixels[offset++] & 0xff;
                    int r = pixels[offset++] & 0xff;
                    data[i] = luminance(r, g, b);
                }
            }*/
            /*else
            {
                throw new ArgumentException("Unsupported image type: " + type);
            }*/
        }

        /// <summary>
        /// Normalizes the contrast.
        /// </summary>
        private void NormalizeContrast()
        {
            int[] histogram = new int[256];
            for (int i = 0; i < m_aData.Length; i++)
            {
                histogram[m_aData[i]]++;
            }
            int[] remap = new int[256];
            int sum = 0;
            int j = 0;
            for (int i = 0; i < histogram.Length; i++)
            {
                sum += histogram[i];
                int target = sum * 255 / m_nPicsize;
                for (int k = j + 1; k <= target; k++)
                {
                    remap[k] = i;
                }
                j = target;
            }

            for (int i = 0; i < m_aData.Length; i++)
            {
                m_aData[i] = remap[m_aData[i]];
            }
        }

        /// <summary>
        /// Writes the edges.
        /// </summary>
        /// <param name="pixels">The pixels.</param>
        private void WriteEdges(int[] pixels)
        {
            //NOTE: There is currently no mechanism for obtaining the edge data
            //in any other format other than an INT_ARGB type BufferedImage.
            //This may be easily remedied by providing alternative accessors.

            /* old Code: edgesImage.getWritableTile(0, 0).setDataElements(0, 0, width, height, pixels);*/

            int i = 0;
            for (int y = 0; y < m_nHeight; y++)
            {
                for (int x = 0; x < m_nWidth; x++)
                {
                    m_SourceImage.SetPixel(x, y, Color.FromArgb(pixels[i]));
                    i++;
                }
            }
        }

        /// <summary>
        /// create a deep copy
        /// </summary>
        /// <returns></returns>
        public override object Clone()
        {
            return new CannyEdgeDetector(this);
        }
    }
}
